struct ECSB2FPTAS <: ECSB2OptimisationAlgorithm end


function optimise_linear_sqrtlinear(instance::CombinatorialInstance, ::ECSB2FPTAS, 
                                    linear::Dict{T, Float64}, sqrtlinear::Dict{T, Float64}, epsilon::Float64, d::Int, verbose::Int; with_trace::Bool=false) where T
  # Solve the problem in polynomial time. Major idea: as the objective function is like w^T x + sqrt(s2^T x), with all x being binaries,
  # can compute all points that maximise the various trade-offs between the two parts of the objective (minus the sqrt part),
  # with a Lagrangian dual-like technique. A maximum value for the dual variable can be derived, and a minimum distance between
  # two trade-offs can be computed. Based on this, use a bisection algorithm to find all the trade-offs. 
  # TODO: Details in the paper (TBD).  

  if with_trace
    runDetails = ECSB2Details()
  end

  # Bound on the number of breakpoints (only useful for the user). solve_linear() maximises the total weight! 
  if verbose >= 1
    x_min = round.(Int, solve_linear(instance, [k => -v for (k, v) in sqrtlinear])) 
    S_min = sum(sqrtlinear[k] * v for (k, v) in x_min)
    x_max = round.(Int, solve_linear(instance, sqrtlinear)) 
    S_max = sum(sqrtlinear[k] * v for (k, v) in x_max)
    K = ceil(Int, (S_max - S_min) / epsilon) 
    
    max_iter = ceil(Int, (d / epsilon) * log2(d / epsilon))
    
    println("About to start iterating. Will find at most $K breakpoints in at most $max_iter iterations.") 
  end

  # Define the discretisation of the lambdas. 
  lambda_incr = (epsilon^2) / (d^2)
  lambda_max = d / epsilon

  # Store the best solution found so far for the nonlinear problem and the set of solutions (helps cut down the search). 
  bestLO = 0.0
  bestNLO = 0.0
  bestLambda = 0.0
  bestSolution = falses(0)
  solutions = Dict{Float64, Vector{T}}()

  # Solve an instance of the linear program and update the best solution so far. 
  function solve_instance(lambda::Float64)
    sol = solve_linear(instance, Dict(arm => linear[arm] + lambda * sqrtlinear[arm] for arm in keys(linear)))
    sum_linear = sum(linear[arm] for arm in sol)
    sum_sqrtlinear = sum(sqrtlinear[arm] for arm in sol)

    solutions[lambda] = sol
    if sum_linear + sqrt(sum_sqrtlinear) > bestNLO
      bestLambda = lambda
      bestSolution = solutions[lambda]
      bestNLO = sum_linear + sqrt(sum_sqrtlinear)
      bestLO = sum_linear + lambda * sum_sqrtlinear
    end
  end

  # The slope is given by the square-root term. 
  slope(lambda::Float64) = sum(sqrtlinear[arm] for arm in solutions[lambda])
  same_slope(lambda1::Float64, lambda2::Float64) = slope(lambda1) == slope(lambda2)

  # Initialise the queue of intervals that must still be explored: the first one is [0, lambda_max]. 
  intervalsToDo = Queue(Pair{Float64, Float64})

  solve_instance(0.)
  solve_instance(lambda_max)

  if ! same_slope(0.0, lambda_max)
    enqueue!(intervalsToDo, 0. => lambda_max)
  end

  # Start iterating over the lambdas and compute the best solution at the same time (initialise with the bounds of lambda). 
  i = 0
  nBreakpoints = 0

  while ! isempty(intervalsToDo)
    interval = dequeue!(intervalsToDo)

    if abs(interval[1] - interval[2]) <= lambda_incr # Very close lambdas? Found a breakpoint! 
      nBreakpoints += 1
    else # Otherwise, continue iterating. 
      lambda = (interval[1] + interval[2]) / 2.
      t0 = now()
      solve_instance(lambda)
      t1 = now()

      if with_trace
        push!(runDetails.solverTimes, (t1 - t0).value)
      end

      # Both ends of the interval have the same slope? Nothing to find here. Otherwise, add the segment. 
      if ! same_slope(interval[1], lambda)
        enqueue!(intervalsToDo, interval[1] => lambda)
      end
      if ! same_slope(lambda, interval[1])
        enqueue!(intervalsToDo, lambda => interval[2])
      end
    end

    i += 1
    if verbose >= 2 && i % 10_000 == 0
      println("Iteration #$i out of $max_iter; total since beginning of loop: $(now() - t0)")
    end
  end

  if verbose >= 1
    println("Done with looping in $i iterations. Took $(now() - t0). Found $nBreakpoints breakpoint$(nBreakpoints == 1 ? "" : "s"). Optimal lambda: $bestLambda. Best objective function: $bestNLO.")
  end

  # Done! 
  if with_trace
    runDetails.nIterations = i
    runDetails.nBreakpoints = nBreakpoints
    runDetails.bestLambda = bestLambda
    runDetails.bestLinearObjective = bestLO
    runDetails.bestNonlinearObjective = bestNLO

    return bestSolution, runDetails
  else
    return bestSolution
  end
end